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Quantum simulation methods based on density-functional theory are currently deemed unfit to 
cope with atomic heat transport within the Green-Kubo formalism, because quantum-mechanical 
energy densities and currents are inherently ill-defined at the atomic scale. We show that, while 
this difficulty would also affect classical simulations, thermal conductivity is indeed insensitive to 
such ill-definedness by virtue of a sort of gauge invariance resulting from energy extensivity and 
conservation. Based on these findings, we derive an expression for the adiabatic energy flux from 
density-functional theory, which allows heat transport to be simulated using ab-initio equilibrium 
molecular dynamics. Our methodology is demonstrated by comparing its predictions with those of 
classical equilibrium and ab-initio non-equilibrium (Miiller-Plathe) simulations for a liquid-Argon 
model, and finally applied to heavy water at ambient conditions. 


Introduction 

Understanding heat transport is key in many fields of 
science and technology, such as materials and planetary 
sciences, energy saving, heat dissipation and shielding, 
or thermoelectric conversion, to name but a few. Heat 
transport in insulators is determined by the dynamics 
of atoms, the electrons following adiabatically in their 
ground state. Simulating atomic heat transport usually 
relies on Boltzmann’s kinetic approach ^ , or on molecu¬ 
lar dynamics (MD), both in its equilibrium (Green-Kubo, 
GK |2H5]) and non-equilibrium flavors. The Boltz¬ 
mann equation only applies to crystalline solids well be¬ 
low melting, whereas classical MD (GMD) bears on those 
materials and conditions that can be modeled by inter¬ 
atomic potentials. Equilibrium ab-initio (AI) MD [3 |S] 
is set to overcome these limitations, but it is still surpris¬ 
ingly thought to be unfit to cope with thermal transport 
because in first-principles calculations it is impossible to 
uniquely decompose the total energy into individual con¬ 
tributions from each atom (excerpted from Ref. [21 • Such 
a unique decomposition is not possible in classical me¬ 
chanics either, because the potential energy of a system 
of interacting atoms can be partitioned into local contri¬ 
butions in an infinite number of equivalent ways. The 
quantum mechanical energy density is also affected by 
a similar indeterminacy. Notwithstanding, the expres¬ 
sion for the heat conductivity derived from any sensible 
energy partitioning or density should obviously be well 
defined, as any measurable quantity must. 

In this work we first demonstrate that the thermal con¬ 
ductivity resulting from the GK relation is unaffected by 
the indeterminacy of the microscopic energy density; we 
then introduce a form of energy density, and a corre¬ 
sponding adiabatic energy flux, from which heat trans¬ 
port coefficients can be computed within the GK for¬ 
malism, using density-functional theory (DFT). Our ap¬ 
proach is validated by comparing the results of equi¬ 
librium AIMD with those of non-equilibrium (Muller- 
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Plathe, MP [2) AIMD and equilibrium GMD simula¬ 
tions for a liquid-Argon model, for which accurate inter¬ 
atomic potentials are derived by matching the forces gen¬ 
erated by them with quantum-mechanical forces com¬ 
puted along the AIMD trajectories. The case of molecu¬ 
lar fluids is finally addressed, and illustrated in the case 
of water at ambient conditions. 

Theory 

According to the GK theory ma, the atomic thermal 
conductivity of an isotropic system is given by: 

^ WkeT^ Jo 

where brackets (•) indicate canonical averages, ks is the 
Boltzmann constant, V and T the system volume and 
temperature, 3q{t) = f(je(r,t) -h (p -h (e)) v(r,t))dr is 
the macroscopic heat flux, jg, v, p, and (e) being the 
energy-current density, local velocity field, and equilib¬ 
rium values of pressure and energy density, respectively 
[nin]. For further reference, we define as diffusive a 
flux that results in a non-vanishing GK conductivity, ac¬ 
cording to Eq. 0 - The integral of the velocity field is 
non diffusive in solids and can be assumed to vanish in 
one-component fluids, because of momentum conserva¬ 
tion. In these cases, as well as in molecular fluids as we 
will see, we can therefore assume that heat and energy 
fluxes coincide. 

Energy is extensive: it can thus be expressed as the 
integral of a density, which is defined up to the divergence 
of a bounded vector field: two densities that differ by such 
a divergence, e(r) and e'(r) = e(r) -|- d ■ p(r), are indeed 
equivalent, in that their integrals differ by a surface term, 
which is irrelevant in the thermodinamic limit, and can 
thus be thought of as different gauges of a same scalar 
field. Energy is also conserved: therefore, for any given 
gauge of its density, e, a corresponding current density, 
je, can be defined so as to satisfy the continuity equation: 

e(r,t)je(r,f) = 0. (2) 

According to Eq. 0 the macroscopic fluxes in two 
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different energy gauges differ by a total time deriva¬ 
tive, which is non-diffusive: Je(t) = Je(t) + P(f), where 
P(<) = f p(r,t)dr. The equality of the corresponding 
heat conductivities results from the following 

Lemma —Let Ji and J2 be two macroscopic fluxes de¬ 
fined for a same system, and 3 12 — Ji -I-J2 their sum. 
The corresponding GK conductivities, Ki, K2, and K12 
satisfy the relation: |ki2 — ki — K2I ^ 2y^KiK2- 

Proof —Let the energy displacement associated with 
the flux Ji be defined as: Di(t) = j,-2 Jo ■ 

The standard Einstein relation HU states that: Ki = 
limi_ioo (|Di(f)P) /L, it follows that: K12 = ki -|- K2 -I- 
limi_,,oo 2 (Di(f) • D2(t)) /t. Canonical averages of prod¬ 
ucts of phase-space functions can be seen as scalar prod¬ 
ucts: the lemma then follows from the Cauchy-Schwartz 
inequality, as applied to the last relation. □ 

The application of the GK methodology to multi- 
component fluids requires some generalizations because 
the presence of multiple atomic species and the exis¬ 
tence of additional hydrodynamical modes (one con¬ 
served number per atomic species) do not permit to iden¬ 
tify the velocity held with the mass-current density, its 
integral with the total momentum, and the heat flux with 
the energy flux. In molecular fluids, however, this iden¬ 
tification can still be done because the integral of the ve¬ 
locity held, while not a constant, is a non-diffusive flux, 
thus not contributing to the heat conductivity. In or¬ 
der to demonstrate this statement, we first define the 
fluxes 3 ab = A ~ where A and B indi¬ 

cate any two atomic species, baIub their stoichiomet¬ 
ric ratio, and Vg = veloci¬ 

ties of all the atoms of a same species S. The integral 
Jq ^AB{t')dt' is equal to the sum of the variations of all 
the AB relative positions in a same molecule, which is 
bound by the sum of the variations of all the AB dis¬ 
tances. Jab is therefore a non-diffusive flux. We have 
N{N — l )/2 such non-diffusive fluxes, N being the num¬ 
ber of species, of which only — 1 are linearly indepen¬ 


dent; furthermore the flux Jm = {Ms is the 

mass of the S-th atomic species) is the total momentum, 
and is thus non-diffusive. We have therefore N inde¬ 
pendent linear combinations of the V5 fluxes that are 
non-diffusive. We conclude that all of them, as well as 
their sum, V(t) = J 2 s'^s{^) ~ V f are also 

non-diffusive. 

In order to derive an expression for the macroscopic 
energy flux appearing in the GK formula, Eq. Q, we 
first multiply the continuity equation, Eq. (H. by r and 
integrate by parts, to obtain the first moment of the time 
derivative of the energy density: 

Je{t) = J e{r,t)rdr. ( 3 ) 

In periodic boundary conditions (PBG) Eq. ([^ is ill- 
defined for the very same reason why macroscopic polar¬ 
ization in dielectrics is so [ 15 ] . In GMD the usual expres¬ 
sion for the energy flux in terms of atomic energies and 
forces [3 is recovered from Eq. ^ by the somewhat 
arbitrary definition: e(r,t) = e/(R, V)( 5 (r — R/), 

where e/ = jM/K/ -h 5 J 2 j^i ^ R-/l)> R- = {R/}. 

and V = {V/} are the atomic energies, positions, 
and velocities, and by reducing the resulting expres¬ 
sion to a boundary-insensitive form. In DFT an en¬ 
ergy density can be defined, which is however inher¬ 
ently ill-determined because of the non-uniqueness of the 
quantum-mechanical kinetic and classical electrostatic 
energy densities nniini. Our previous analysis demon¬ 
strates that, in spite of previous worries to the contrary, 
the transport coefficients derived from a DFT energy 
density through the GK formula, Eq. 0 , are well de¬ 
fined, provided a macroscopic energy flux can be com¬ 
puted from Eq. ^ in PBG. To this end, among many 
equivalent gauges, we choose to represent the DFT total 
energy as the integral of the density: 


eoFrir) = '^ 6 {r -Ri)e°j + Re'^ip*{r){HKS<Pv{r)) - ^p{r)vH{r) + {exc{r) - vxcir)) p{r), ( 4 ) 

I V 


where ej = ^MjVj + wj are bare ionic energies; Mj, 
Zj, and ui/ = I being ionic masses, 

charges, and electrostatic energies, respectively; the elec¬ 
tron charge is assumed to be one; Hks is the instanta¬ 
neous Kohn-Sham (KS) Hamiltonian, (p„’s its occupied 
eigenfunctions, and p(r) = |(/?„(r)p the ground-state 

electron-density distribution; vb and vxc Hartree 
and exchange-correlation (XG) potentials, and exc is a 
local XG energy per particle, defined by the relation: 


Exc = / exc[p](i’)p(i*)c^i'n The energy density of Eq. 
Q depends on time through atomic positions and veloc¬ 
ities and KS orbitals. Inserting its time derivative into 
Eq. (|^ and using the Born-Oppenheimer (BO) equa¬ 
tions of motion for the nuclei {Mj'Vj = —dEBPr/dRi), 
the resulting adiabatic energy flux can be expressed as: 

Je = JifS + Jh + Jq + Jo + Jxc- (5) 


^ ^XC is also to some extent ill-defined, in that any XC densities 
resulting in a same integral should be considered as equivalent. 
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The five fiuxes in Eq. ([^ are defined as: 

^KS = X! (^{^v\rHKs\^v) + £v{^v\A^v)^ , (6) 

V 

^ j VH{r)yvH{r)dr, (7) 

vj 

Jc^eM + E (Rz-RL)(Vi.Viraj)], (9) 

/ L^I 

T =/° noi 

[-/p(r)p(r)(9eGGA(r)rfr (GGA), 

where e„ in Eq. is the w-th eigenvalue of the KS 

Hamiltonian; V and V/ = in Eqs. (Tjjfil in- 


' or ^ ot^T ^ A H 

dicate the gradients with respect to the argumenTra the 
function and to the /-th atomic position, respectively; 
the simbol vq in Eq. Q indicates the (possibly non¬ 
local) ionic (pseudo-) potential acting on the electrons; 
finally, “LDA” and “GGA” in Eq. (101 indicate the local- 


density m and generalized-gradient |18) approximations 
to the XG energy functional, respectively, and decGA the 
derivative of the GGA XG local energy per particle with 
respect to density gradients. Eq. (§ can be derived from 
Eqs. (IB with some tedious but straightforward algebra 
(see Methods). The last four terms on its right-hand side, 
Eqs. (|7 lOl, are manifestly boundary-insensitive, whilst 
the first, Eq. ([fi]), is not, because the position operator 
appearing therein is ill-defined in PEG. Within the adia¬ 
batic time evolution that is assumed in AIMD, however, 
the time derivative of a KS orbital, as well as its product 
with the KS Hamiltonian, are orthogonal to the orbital 
itself in the “parallel transport” gauge where KS orbitals 
are real |l9j[20]|^ {'^v\'^v) = 0 cOid l^y\HKsWv) = 0. 
Therefore, in order to evaluate Eq. ([6|, one only needs 
the projection of r|(/?„) onto the manifold orthogonal to 
which is well defined in PEG. Actually, by expanding 
ify in the basis of the eigenstates of the instantaneous 
KS Hamiltonian m, one sees that only the projection of 
r|(^„) onto the empty-state manifold, \(p^) = PcX°'\(py), 
contributes to J_ffs, where Pc = 1—is 
the Qf-th Gartesian component of r. Using the standard 
prescription adopted in density-functional perturbation 
theory (DEPT), such a projection can be computed by 
solving the linear equation m- 


{HkS - £v)\^v) = Pc[HKS,X°‘]\ipPl, 


( 11 ) 


where the ill-definedness of the solution, due to the sin¬ 
gularity of the left-hand side, is lifted by enforcing its 


^ The concept of gauge for the quantum-mechanical representation 
of molecular orbitals should not be confused with that introduced 
in this paper for the energy density. 


orthogonality to the occupied-state manifold. In terms 
of the (py’s Eq. (U reads: 


JPs = E + £a(^a|^“)) ■ (12) 

V 


The flux in Eq. (121 is not manifestly invariant with 
respect to the arbitrary choice of the zero of the one- 
electron energy levels. A shift of the energy zero by 
a quantity Ae results in a shift of the Kohn-Sham en¬ 
ergy flux: Jfs ^ Jfs -h {{(Pv\^v) + {^v\(Pv)) = 

J'^g Ae , where Jp is the adiabatic electronic macro¬ 
scopic flux introduced in Ref. m- The electronic cur¬ 
rent is the difference between the total charge current 
and its ionic component: the first is by definition non- 
diffusive in insulators, while the second is so in mono- 
atomic and molecular systems, as we have seen when 
discussing the latter. We conclude that the electronic 
flux is non-diffusive in insulators, thus not contributing 
to their heat conductivity and lifting the apparent inde¬ 
terminacy of Eq. (121. 


Numerical simulation 

The methodology presented above has been implemented 
in the Quantum ESPRESSO suite of computer codes 
|22| : a Gar-Parrinello (GP) j7j AIMD trajectory is first 
generated using the cp. x code; the energy flux is then 
evaluated along this trajectory according to Eqs. ffl 
by an add-on to the pw. x code implemented using sev¬ 
eral DEPT routines borrowed from the ph.x code; the 
thermal conductivity is finally computed from the GK 
relation, Eq. ([^, or the equivalent Einstein relation la¬ 
in order to demonstrate this methodology, we com¬ 
pare its predictions with those from GMD j23] for a sys¬ 
tem whose DFT EO energy surface can be accurately 
mimicked by pair potentials. Not aiming at a real¬ 
istic description of any specific system, but rather at 
the ease and accuracy of the classical representation of 
the DFT EO surface, we choose liquid Argon and use 
the LDA XG functional, in spite of the well known in¬ 
ability of the latter to capture dispersion forces. This 
reference system will be dubbed “LDA-Ar”. KS or¬ 
bitals are treated within the plane-wave (PW) pseudo¬ 
potential (PP) method [51]. Our model consists of 108 
atoms in a periodically repeated cubic supercell with 
an edge of 33 a.u., corresponding to a density of 1.34 
g/cnP. AIMD trajectories were generated via the Gar- 
Parrinello dynamics |7| for 100 picoseconds (ps), using a 
time step of 0.242 femtoseconds (fs) and a fictitious elec¬ 
tronic mass of 1000 electronic masses, at two different 
temperatures, T = 250 and 400 K. The fictitious elec¬ 
tronic temperature was monitored and checked not to 
be subject to any significant drift. The EO energy sur¬ 
face was modeled with a sum of classical pair potentials 
of the form V{r) = P 2 {r)e~°‘^, where P 2 is a second- 
order polynomial, whose parameters were determined in¬ 
dependently for each temperature by a least-square fit 
of the classical vs. quantum-mechanical forces computed 
along the AIMD trajectory. Self-diffusion coefficients of 
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(10.8±0.1), and (15.6±0.2) x 10“®cm^/s were estimated 
along the two AIMD trajectories, in close agreement with 
the CMD values (10.3±0.1), and (15.8±0.2) x 10“®cm^/s, 
thus confirming the quality of the classical model. Radial 
distribution functions computed from AIMD and CMD 
trajectories were also found to be very similar. 




FIG. 1. Time correlations of the energy flux in LDA Ar. 
(a) C4t) = . J,(0))[10i®mWK-is-']. (b) 

K{t) = fg Ce{t')dt' [mWm~^K“^] (see text). Blue: AIMD 
(100 ps). Orange, CMD (100 ps). Green, CMD (1000 ps). 
The shaded areas depict statistical errors as estimated from 
a block analysis of our MD trajectories. 

In Fig. we compare the time correlation func¬ 
tions of the energy flux in LDA Ar, as computed from 
AIMD and CMD at T = 250 K. The CMD and AIMD 
correlation functions differ not quite because they cor¬ 
respond to different systems—which are actually close 
enough as to have very similar equilibrium and diffusion 
properties—as because the AIMD and CMD fluxes de¬ 
rive from a different unpacking of the total energy into 
local contributions. In Fig. [^3 we display the integrals 
= So 3 Vknr^ • Je(0))dt'; the AIMD and CMD 
heat conductivities, k = limt_,,oo ^{t), coincide within 
statistical errors with each other and with the CMD value 
evaluated from a 1-ns-long simulation: (103 ± 5, 100 ± 6, 
and 104 ± 2) [mW K“^m“^], respectively. A similar 
level of agreement is obtained for the other temperature, 
T = 400K (118±8, 112±7, and 110±2) [mW 

In order to further validate these results, we have re¬ 
computed the thermal conductivities of our LDA-Argon 
model, using non-equilibrium (MP) AIMD |5]. A detailed 
comparison of GK vs. MP AIMD for heat-transport sim¬ 


ulations is out of the scope of the present paper, and 
we have limited ourselves to two MP simulations, aimed 
at mimicking the physical conditions of the GK AIMD 
simulations reported above, and performed using mini¬ 
mal simulation settings: we used (2x2x5) supercells, 
where the notation indicates multiples of a 4-atom cu¬ 
bic unit cell, thus resulting in 80-atom tetragonal super¬ 
cells whose size was chosen so as to result in the same 
mass density of 1.34 g/cm^ as used before. MP simu¬ 
lations were performed by subdividing the supercell in 
eight equally spaced layers stacked along the c axis and 
by swapping the velocities of the hottest atom in the 
cool region and the coolest atom in the hot region ev¬ 
ery picosecond. Rather long simulations (> 360 ps) were 
necessary to achieve an acceptable statistical accuracy, 
resulting in estimated thermal conductivities of 94 ± 13 
and 109 ± 11 [mW at the temperatures of 287 

and 423 K, respectively. Our GK and MP AIMD re¬ 
sults are compared in Fig. witnessing to a convincing 
validation of our approach based on the Green-Kubo for¬ 
malism. 



FIG. 2. Gomparison of the heat conductivities of our LDA- 
Ar model, as estimated from Green-Kubo and Miiller-Plathe 
ab-initio molecular dynamics. Units are [mW K“^m~^]. Or¬ 
ange: equilibrium (GK) molecular dynamics; the two dots 
indicate the estimates from our simulations, the straight line 
their linear inter-/extrapolation. Statistical errors, as esti¬ 
mated by a block analysis of our MD trajectories, are indi¬ 
cated by error bars or by shaded areas, where relevant. 

We have applied our newly developed method to com¬ 
pute the heat conductiviy of heavy water at ambient con¬ 
ditions. We have generated a 90-ps long AIMD trajec¬ 
tory for a system of 64 heavy-water molecules in a cubic 
supercell with an edge of 23.46 a.u., corresponding to 
the experimental density of 1.11 g/cm^, and at an es¬ 
timated temperature T = 385 K, using the PBE XC en¬ 
ergy functional [18] and the PW-PP method as above 
|24| . A time step of 0.0726 fs and a fictitious elec¬ 
tronic mass of 340 electron masses were used in this case. 
The resulting self-diffusion coefficient was estimated to 
(2.6 ± 0.2) X 10“® cm^/s, to be compared with an exper¬ 
imental value of 2.0 x 10“®cm^/s at T = 298 K, follow¬ 
ing the common practice of comparing experimental data 
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FIG. 3. Energy displacement in liquid heavy water at ambient 
conditions. D^{t) [10~^^mJm“^K“^] was evaluated from the 
modified flux J* = Je + A*V at T = 385 K with a GGA 
XG functional (see text); the dashed line indicates a linear 
fit to the large-time behavior of the curve. Inset: integral 
in the GK relation, K{t) Eq. Q, as a function of the upper 
limit of integration (see caption to Fig. 1); the dashed line 
indicates the value of the thermal conductivity obtained from 
the Einstein relation, i.e. from the slope of above linear fit. 
The shaded areas depict statistical errors, as estimated from 
a block analysis of our MD trajectories. 


, as a function of t in 


for water at ambient conditions with AIMD-PBE simula¬ 
tions performed at ^ 400 K . The power spectrum of 
the computed energy flux is characterized by three rela¬ 
tively narrow peaks in correspondence to the intramolec¬ 
ular vibrational modes |26| . resulting in long-lived high- 
frequency oscillations in the integrand of Eq. (Q, that 
plague the evaluation of the integral as a function of the 
upper limit of integration well beyond the time where the 
noise of the integrand becomes larger than the amplitude 
of its oscillations. As the computation of transport co¬ 
efficients from the Einstein relation m is less affected 
by the high-frequency components of the power spectrum 
m, this ailment is alleviated by evaluating the heat con¬ 
ductivity as the slope of the energy squared displaeement, 

^ei^) = eVksT^ fo Je(i')dt' 

the large-time limit. A direct application of this tech¬ 
nique is however not possible as the long-time behavior 
of the energy squared displacement does not allow us to 
extrapolate a straight line before it becomes too noisy to 
analyze. This state of affairs indicates the existence of a 
slowly decaying mode in the energy-ffux correlation func¬ 
tion, possibly correlated with a non-diffusive flux. As we 
have seen, the total velocity V is such a non-diffusive 
flux. The value of the corresponding GK conductivity, 
Eq. ([^, however, goes to zero very slowly as a func¬ 
tion of the upper limit of integration. This suggests that 
the slow convergence of the heat conductivity of water as 
estimated from the slope of the energy squared displaca- 
ment as a function of time, is possibly due to large cor¬ 
relations existing between the energy flux and the total 
velocity. We have therefore decided to analyze, instead 
of Je, the modified flux J* = Je + A*V, where A* has 


been fixed in such a way as to minimize the correlations 
between J* and V. Fig. 2 displays the squared energy 
displacement computed from J* as a function of time and 
demonstrates that a constant slope can indeed be identi¬ 
fied in the long-time limit, giving a value for the heat con¬ 
ductivity of heavy water of 740 ± 120 mWm“^K“^, to be 
compared with an experimental value of 606 mWm“^K“^ 
and 595mWm“^K“^ for light and heavy water respec¬ 
tively at ambient conditions ESI- The inset displays the 
behavior of K(t) (see caption to Fig. 1) as a function of 
the upper limit of integration in the GK formula, indi¬ 
cating that a direct use of Eq. (Q would be extremely 
difficult in this case. A more detailed error analysis and 
a systematic extension of this study to different isotopic 
compositions and other conditions of temperature and 
pressure is currently in the works. 

Conclusions 

We believe that the discussion presented in this work will 
elucidate the scope of a number of assumptions that, al¬ 
though routinely made in the classical simulation of heat 
transport, have never been fully clarified, thus hamper¬ 
ing their generalization to quantum simulations. We are 
confident that the resulting new methodology will have 
an impact on important problems where other methods 
may fail, such as e.g. liquids and glasses, particularly at 
extreme conditions of temperature and pressure. 

Methods 

In order to derive Eqs. ([5 
which we rewrite as: 

eoFrir) = e/fs(r) + eo(r) -f e//(r) -b exci^), (13) 
where 

eifs(r) = Re Y,<P*v{r){HKsM^)), (14) 

V 

eo(i’) =^^(r-R/) , (15) 

eff(r) =-^p(r)ui/(r), and (16) 

exc{r) = {exc{r)-vxc{r)) p{r), (17) 

exc is a local XG energy per particle, defined by the 
relation 


10), we start from Eq. 0. 


Exc 


= j exc[p](r)p(i-)dr, 


(18) 


and the XG potential vxc is 


i^xc(r) 


SExc 

5p{v) 


= exc(r) 


5exc{E) 

Sp{r) 


p{r')dr'. 


(19) 


In the LDA, exc is a function of the local density, 
whereas in the GGA it is a function of the local density 
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and density gradients: 


= f^LDA{pir)), 

<^xc^[pK^) = eGGA(p(r), Vp(r)). 


We now proceed to computing the first moments of the 
time derivatives of the above four densities, according to 

(20) III order to simplify the notation, the time de¬ 
pendence of the various quantities will be omitted. Let’s 

( 21 ) 

start with the Kohn-Sham energy density, Eq. (1T3. 


eKs(r) = + >pl{r)HKS‘Pv{r) + <pl{r)HKS‘Pv{r)'^ 

V 

= eKsir:) + eo(r) + e'ij(r) -h W, 


( 22 ) 

(23) 


where 


SKsir) = [e„(/5*(r)(/5„(r) -h ipl{r)H ks^P vir) , (24) 

V 

e'oir) = Y ‘P*vi^)^o<Pnir), (25) 


eff(i’) ='i’ff(i’)p(r), and 
= vxci'r)p{r). 


(26) 

(27) 


The macrosopic flux deriving from eKS^ Eq- (241, is the 


“Kohn-Sham” flux of Eq. 

J reKs{^)dr = Jks- (28) 


The other three terms, Eqs. (25 271 result from the 


external-, Hartree-, and XC-potential contributions to 
the time derivative of the KS Hamiltonian (third term 
in Eq. 22). The corresponding fluxes combine with the 


fluxes originating from the energy densities of Eqs. (15 
17), as explained below. 


The first moment of the “ionic potential” energy- 
density derivative, Eq. (25), reads: 


reo(r)dr = ^ (</?.„ |rt}ofy„) 

V 

= Y I"* ■ ^I^o)\‘Pv) 

vj 


= ■ '^^'^o)| ^n) + R/(v?-«|(V/ • V/Uo)fy„) 

V,I 

= J'o-Y^i (V/-FfO, 


(29) 


where Jq is the flux of Eq. ([^, and Ff is the electronic from the first moment of the “ionic” energy density, Eq. 
(Hellmann-Feynman) contribution to the force acting on (151. 
the J-th atom. The corresponding (second) term in the 

energy flux of Eq. (p^ is ill-defined in PBC but, as we The time derivative of the “ionic” energy density, Eq. 


will see shortly, it cancels with a similar term coming (15), reads: 


eo(r) = 


I 


Vi ■ ViS{r - R/) + 6{r - Rj) ( MjVi • V/ + ^ Vj • Vjwi 

j^i 


(30) 


We now use Newton’s equations of motion {MiVi = F/, 
where F/ is the force acting on the I-th atom), and split 


F/ into an electronic (Hellmann-Feynman) contribution, 
plus a sum of pair-wise electrostatic terms, F/ = Ff — 
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^ V/Wj, to obtain: 


reo 


3o(r)dr = ^ [e?V/ + R, (f, • V, 4 
/ 

= E[e?V/ + R/ (Ff-Yi) 

I 

= E[e?V/ + R/ (Ff-V,) 
/ 

= Jo 4-ER/ (Ff-V,), 


.747 

+ R-7 E (^''' ■ ~ IWJ 

J^I 

+ ^(R7-R,7)(Vj.V,777;7)] 

J^I 


(31) 


r 


where Jq is the energy flux of Eq. ^ and the third step 
follows from the second by interchanging the dummy in- 
deces of one of the two sums over I and J. As antici¬ 
pated before, the second term on the right-hand side of 
Eq. (311, which is ill-defined in PBC, cancels a similar 


term in Eq. (291, leaving all the surviving terms well 


defined. We summarize Eqs. (291 and (311 as: 


y r[eo(r) + eo(r)]dr = Jq 4-Jq, 


(32) 


where Jq and Jq are the energy fluxes of Eqs. (§ and 
respectively. 


We then combine the time derivative of the “Hartree” 


energy density, Eq. (161, with the “Hartree-potential” 


energy-density derivative, Eq. (26l: 


677 (r) = 677 (r) + (33) 

= ^(«77(r)p(r) -p(r)u77(r)) 

= - VH{r)AvH{r)) 

= ■ {vH{r)\7vH{r) - VH{r)YvH{r)) (34) 


Multiplying Eq. (341 by r and integrating by parts, one 
obtains: 


J 77 = y re 77 (r)dr 

^ ^ / '^H{r)VvH{r)dr, 


(35) 


which is Eq. 0- 

We Anally address the first moments of the time deriva¬ 
tive of the “XC” energy density, Eq. ( [T7| , and of the 
“XC-potential” energy-density derivative, Eq. ( [2^ . We 
define: 


exc(r) = exc{r) + excir) 

= (exc(r) - 'yxc(i’))p(r) -4 excir)pir) 


Sexci^ ) , I 

c f. )dr , 

dp{r) 


(36) 


r 


which derives from the definition of the XC potential. The first moment of Eq. (36) reads: 
Eq. (191, and from the chain rule as applied to the time 


derivative of exc' 


= J • 


(37) 


Jxc = J rexc{'r)dr 

= jiy- F)p(r)p(r0 ^^yyj'^ drdr'. 


(38) 





















In the LDA, because of the local dependence of exc on 


the electron density, the functional derivative in Eq. (38) 


is proportional to ^(r —r'), thus making the integral van¬ 
ish. In the GGA Eq. ( [^ gives: 


Sp{r') 


= eGG.4(i')<5(r- !■') 

+ daeGGA{r)Va6{r - r'), (39) 


where 


'GGA 


(r) = 


^ OeaaA(pyp) ^ d^ecGAir) = 

p=p(r) 

. The first term on the right-hand side 

P=p(r) 

of Eq. ( |39[ ) does not contribute to the XG energy flux as 
in the LDA. By inserting the second term into Eq. (381, 


duGGAipyp) 

dVcP 


one finally arrives at the expression for the XG energy 
flux of Eq. (101, thus completing the derivation of Eqs. 

This rather unwieldy, but all in all straightfoward, 
derivation is visually summarized in Fig. 


r 

^0 


^KS 

u 


Cq ^KS ^xc^h 


Jo “1“ Jq ^ KS 


I 


'H ^XC 


1 




L 


'H ^XC 


^xc 

J 


FIG. 4. Conceptual flow of the derivation of the various 
components of the macroscopic energy flux, Eqs. ( |6|10[), fro m 
the definition of a microscopic energy density, Eqs. ( 13|17 1 
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